---
title: "Empirical evidence of mental health risks posed by climate change"
subtitle: Main text replication code
author: Nick Obradovich, Robyn Migliorini, Martin P. Paulus, and Iyad Rahwan
output: pdf_document
---

<!--
Rscript -e "rmarkdown::render('./main_text_replication_code.Rmd')"
-->

```{r opts, echo=FALSE, message=FALSE, warning=FALSE, include=FALSE}
knitr::opts_chunk$set(fig.path='.', warning=F, message=F, fig.retina=NULL, cache=T, cache.lazy=F, autodep=T, echo=F, eval=T)
```

```{r packages}
# packages----
library(bit64) # long integers
library(ggthemes) #themes for plotting
library(sandwich) #robust se estimator
library(lmtest) #lm functions
library(ordinal) #ologit funcitons
library(memisc) #import stata data
library(ggmap) #make nice plots with ggplot mapping
library(ggplot2) #nice plots
library(Hmisc) #multipurpose package
library(RColorBrewer) #make nice colors
library(grid) #plotting package
library(MASS) #multipurpose package
library(gplots) #nice plotting
library(lfe) #run linear fixed effects models
library(stargazer) #output nice tables
library(lubridate) #deal with dates
library(zoo) #deal with dates
library(reshape2) #powerful reshape package
library(sp) #general spatial class package
library(gstat) #geostat package
library(spacetime) #store data in proper spatial/temporal class
library(raster) #convert to raster and vice versa
library(maptools) #plot and modify shapefiles
library(rgeos) #for use with maptools for shapefiles
library(rgdal) #for reading in shapefiles
library(RSAGA) #needed for spatial downscaling
library(RCurl) #also needed for spatial downscaling
library(dplyr) #for data frame manipulation functions
library(plyr) #for additional manipulation functions
library(data.table) #for data.table functions
library(foreign) #for reading in .dta files
library(parallel) #for parallel processing functions
library(foreach) #for plyr parallel
library(doMC) #run multicore processes
library(caTools) #functions for fast running mean and sd
library(ncdf4) #tools for netCDF packages
library(stringr) #string tools
library(knitr) #knitting
library(pander) #pandering
library(rasterVis) # plotting rasters
library(GGally) # for correlation plots
library(geosphere) # spherical distance
library(Rcpp) # for C++
library(RcppArmadillo) # for C++ (and Conley errors)
library(matrixStats) # matrix calculations
library(ggExtra) # for marginal histograms
library(gridExtra) # for multiplots
library(viridis) # for colors
library(Cairo) # pdf device
library(gtable) # for gtable filter
library(devtools)
#devtools::install_github("wilkelab/cowplot")
library(cowplot)
library(captioner)

# set cores for parallel processing----
registerDoMC(20)
```

```{r functions}
# functions for grabbing p-values and coefficients from summaries of felm models----
coefs <- function(reg, name) {
    if (round(reg$coefficients[rownames(reg$coefficients)==name, 1], 3) == 0) {
        format(round(reg$coefficients[rownames(reg$coefficients)==name, 1], 4), scientific=FALSE) 
    }
    else round(reg$coefficients[rownames(reg$coefficients)==name, 1], 3)
}
pval <- function(reg, name) {
  if (round(reg$coefficients[rownames(reg$coefficients)==name, 4], 3) >=0.001) 
  {round(reg$coefficients[rownames(reg$coefficients)==name, 4], 3)} 
  else paste("< 0.001")
}

# captioner setting----
fig_num <- captioner(prefix="Fig.")
tab_num <- captioner(prefix="Table")
eq_num <- captioner(prefix="Equation")
sfig_num <- captioner(prefix="Fig. S", auto_space=FALSE)
stab_num <- captioner(prefix="Table S", auto_space=FALSE)
seq_num <- captioner(prefix="Equation S", auto_space=FALSE)

# binned plot----
source('./binned_plot_function.R')
```

```{r data}
load('./climate_mental_health.Rdata')
```

```{r models}
# contemporaneous regression----
fig1 <- summary(felm(any_ment ~ cutprcp + cuttmax + avg.trange + avg.humid + avg.cloud + avg.wind
          | time + mmsa + stateyr
          | 0 
          | state
          , data=brfss,
          exactDOF=TRUE))

# heterogeneous effects regressions----
fig2.poor <- summary(felm(any_ment ~ cutprcp + cuttmax + avg.trange + avg.humid + avg.cloud + avg.wind
            | time + mmsa + stateyr
            | 0
            | state
            , data=brfss[income2 <= 22500,],
            exactDOF=TRUE))
fig2.rich <- summary(felm(any_ment ~ cutprcp + cuttmax + avg.trange + avg.humid + avg.cloud + avg.wind
            | time + mmsa + stateyr
            | 0
            | state
            , data=brfss[income2 >= 90000,],
            exactDOF=TRUE))
fig2.female <- summary(felm(any_ment ~ cutprcp + cuttmax + avg.trange + avg.humid + avg.cloud + avg.wind
            | time + mmsa + stateyr
            | 0
            | state
            , data=brfss[female == 1,],
            exactDOF=TRUE))
fig2.male <- summary(felm(any_ment ~ cutprcp + cuttmax + avg.trange + avg.humid + avg.cloud + avg.wind
            | time + mmsa + stateyr
            | 0
            | state
            , data=brfss[female == 0,],
            exactDOF=TRUE))
fig2.poorfemale <- summary(felm(any_ment ~ cutprcp + cuttmax + avg.trange + avg.humid + avg.cloud + avg.wind
            | time + mmsa + stateyr
            | 0
            | state
            , data=brfss[income2 <= 22500 & female == 1,],
            exactDOF=TRUE))
fig2.richmale <- summary(felm(any_ment ~ cutprcp + cuttmax + avg.trange + avg.humid + avg.cloud + avg.wind
            | time + mmsa + stateyr
            | 0
            | state
            , data=brfss[income2 >= 90000 & female == 0,],
            exactDOF=TRUE))

# long differences regression----
pre.years <- 2002:2006
post.years <- 2007:2011

## respondent pre-period averages
r.pre <- brfss[year(time) %in% pre.years,
               .(any_ment = unlist(lapply(.SD, mean, na.rm=T)),
                 years = uniqueN(year(time)),
                 subs = uniqueN(.I)),
               by=.(mmsa, state, lat, lon),
               .SDcols = c("any_ment")]
setnames(r.pre, c("any_ment", "years", "subs"), paste0(c("any_ment", "years", "subs"), '.pre'))
setkey(r.pre, mmsa)

## city pre-period averages
prism[, trange := tmax - tmin]
c.pre <- prism[year(time) %in% pre.years,
               lapply(.SD, mean, na.rm=T),
               by=.(mmsa),
               .SDcols = c("tmax","trange","prcp","cloud","humid","wind")]
setnames(c.pre, c("tmax","trange","prcp","cloud","humid","wind"), paste0(c("tmax","trange","prcp","cloud","humid","wind"), '.pre'))
setkey(c.pre, mmsa)

## merge pre
pre <- merge(r.pre, c.pre)
rm(r.pre, c.pre)

## respondent post-period averages
r.post <- brfss[year(time) %in% post.years,
                .(any_ment = unlist(lapply(.SD, mean, na.rm=T)),
                  years = uniqueN(year(time)),
                  subs = uniqueN(.I)),
                by=.(mmsa, state, lat, lon),
                .SDcols = c("any_ment")]
setnames(r.post, c("any_ment","years","subs"), paste0(c("any_ment","years","subs"), '.post'))
setkey(r.post, mmsa)

## city post-period averages
c.post <- prism[year(time) %in% post.years,
                lapply(.SD, mean, na.rm=T),
                by=.(mmsa),
                .SDcols = c("tmax","trange","prcp","cloud","humid","wind")]
setnames(c.post, c("tmax","trange","prcp","cloud","humid","wind"), paste0(c("tmax","trange","prcp","cloud","humid","wind"), '.post'))
setkey(c.post, mmsa)

## merge post
post <- merge(r.post, c.post)
rm(r.post, c.post)

## merge
setkey(pre, mmsa, state, lat, lon)
setkey(post, mmsa, state, lat, lon)
diff <- merge(pre, post) # merge together

## calculate total number of years that enter sample
diff[, years := years.pre + years.post]
diff[, subs := subs.pre + subs.post]

## list of variables to difference
v <- gsub('.pre', '', grep('.pre', names(pre), value=TRUE))
invisible(
    foreach::foreach (l=v) %do% {
        diff <- diff[, paste0(l, '.diff') := get(paste0(l, '.post')) - get(paste0(l, '.pre'))]
    } # generate all differences variables
)

diff[, c('years.diff','subs.diff') := NULL]
diff <- diff[!is.nan(tmax.diff),] # remove obs with NaN
diff[, tmax.diff2 := tmax.diff^2]

ld <- summary(felm(any_ment.diff ~ tmax.diff + prcp.diff + trange.diff + cloud.diff + humid.diff + wind.diff | 0 | 0 | state, data=diff, weights=diff[, subs]))
ld.full <- summary(felm(any_ment.diff ~ tmax.diff + prcp.diff + trange.diff + cloud.diff + humid.diff + wind.diff | 0 | 0 | state, data=diff[years==10,], weights=diff[years==10, subs])) # cities with all 10 years in sample
ld.partial <- summary(felm(any_ment.diff ~ tmax.diff + prcp.diff + trange.diff + cloud.diff + humid.diff + wind.diff | 0 | 0 | state, data=diff[years<10,], weights=diff[years<10, subs])) # cities will less than 10 years in sample

# long-diff spring----
pre.years <- 2002:2006
post.years <- 2007:2011

## respondent pre-period averages
r.pre <- brfss[(month(time) %in% c(3,4,5)) & (year(time) %in% pre.years),
               .(any_ment = unlist(lapply(.SD, mean, na.rm=T)),
                 years = uniqueN(year(time)),
                 subs = uniqueN(.I)),
               by=.(mmsa, state, lat, lon),
               .SDcols = c("any_ment")]
setnames(r.pre, c("any_ment", "years", "subs"), paste0(c("any_ment", "years", "subs"), '.pre'))
setkey(r.pre, mmsa)

## city pre-period averages
prism[, trange := tmax - tmin]
c.pre <- prism[(month(time) %in% c(3,4,5)) & (year(time) %in% pre.years),
               lapply(.SD, mean, na.rm=T),
               by=.(mmsa),
               .SDcols = c("tmax","trange","prcp","cloud","humid","wind")]
setnames(c.pre, c("tmax","trange","prcp","cloud","humid","wind"), paste0(c("tmax","trange","prcp","cloud","humid","wind"), '.pre'))
setkey(c.pre, mmsa)

## merge pre
pre <- merge(r.pre, c.pre)
rm(r.pre, c.pre)

## respondent post-period averages
r.post <- brfss[(month(time) %in% c(3,4,5)) & (year(time) %in% post.years),
                .(any_ment = unlist(lapply(.SD, mean, na.rm=T)),
                  years = uniqueN(year(time)),
                  subs = uniqueN(.I)),
                by=.(mmsa, state, lat, lon),
                .SDcols = c("any_ment")]
setnames(r.post, c("any_ment","years","subs"), paste0(c("any_ment","years","subs"), '.post'))
setkey(r.post, mmsa)

## city post-period averages
c.post <- prism[(month(time) %in% c(3,4,5)) & (year(time) %in% post.years),
                lapply(.SD, mean, na.rm=T),
                by=.(mmsa),
                .SDcols = c("tmax","trange","prcp","cloud","humid","wind")]
setnames(c.post, c("tmax","trange","prcp","cloud","humid","wind"), paste0(c("tmax","trange","prcp","cloud","humid","wind"), '.post'))
setkey(c.post, mmsa)

## merge post
post <- merge(r.post, c.post)
rm(r.post, c.post)

## merge
setkey(pre, mmsa, state, lat, lon)
setkey(post, mmsa, state, lat, lon)
diff <- merge(pre, post) # merge together

## calculate total number of years that enter sample
diff[, years := years.pre + years.post]
diff[, subs := subs.pre + subs.post]

## list of variables to difference
v <- gsub('.pre', '', grep('.pre', names(pre), value=TRUE))
invisible(
    foreach::foreach (l=v) %do% {
        diff <- diff[, paste0(l, '.diff') := get(paste0(l, '.post')) - get(paste0(l, '.pre'))]
    } # generate all differences variables
)

diff[, c('years.diff','subs.diff') := NULL]
diff <- diff[!is.nan(tmax.diff),] # remove obs with NaN

ld.spr <- felm(any_ment.diff ~ tmax.diff + prcp.diff + trange.diff + cloud.diff + humid.diff + wind.diff | 0 | 0 | state, data=diff, weights=diff[, subs])
ld.full.spr <- felm(any_ment.diff ~ tmax.diff + prcp.diff + trange.diff + cloud.diff + humid.diff + wind.diff | 0 | 0 | state, data=diff[years==10,], weights=diff[years==10, subs]) # cities with all 10 years in sample
ld.partial.spr <- felm(any_ment.diff ~ tmax.diff + prcp.diff + trange.diff + cloud.diff + humid.diff + wind.diff | 0 | 0 | state, data=diff[years<10,], weights=diff[years<10, subs]) # cities will less than 10 years in sample

# long-diff summer----
pre.years <- 2002:2006
post.years <- 2007:2011

## respondent pre-period averages
r.pre <- brfss[(month(time) %in% 6:8) & (year(time) %in% pre.years),
               .(any_ment = unlist(lapply(.SD, mean, na.rm=T)),
                 years = uniqueN(year(time)),
                 subs = uniqueN(.I)),
               by=.(mmsa, state, lat, lon),
               .SDcols = c("any_ment")]
setnames(r.pre, c("any_ment", "years", "subs"), paste0(c("any_ment", "years", "subs"), '.pre'))
setkey(r.pre, mmsa)

## city pre-period averages
prism[, trange := tmax - tmin]
c.pre <- prism[(month(time) %in% 6:8) & (year(time) %in% pre.years),
               lapply(.SD, mean, na.rm=T),
               by=.(mmsa),
               .SDcols = c("tmax","trange","prcp","cloud","humid","wind")]
setnames(c.pre, c("tmax","trange","prcp","cloud","humid","wind"), paste0(c("tmax","trange","prcp","cloud","humid","wind"), '.pre'))
setkey(c.pre, mmsa)

## merge pre
pre <- merge(r.pre, c.pre)
rm(r.pre, c.pre)

## respondent post-period averages
r.post <- brfss[(month(time) %in% 6:8) & (year(time) %in% post.years),
                .(any_ment = unlist(lapply(.SD, mean, na.rm=T)),
                  years = uniqueN(year(time)),
                  subs = uniqueN(.I)),
                by=.(mmsa, state, lat, lon),
                .SDcols = c("any_ment")]
setnames(r.post, c("any_ment","years","subs"), paste0(c("any_ment","years","subs"), '.post'))
setkey(r.post, mmsa)

## city post-period averages
c.post <- prism[(month(time) %in% 6:8) & (year(time) %in% post.years),
                lapply(.SD, mean, na.rm=T),
                by=.(mmsa),
                .SDcols = c("tmax","trange","prcp","cloud","humid","wind")]
setnames(c.post, c("tmax","trange","prcp","cloud","humid","wind"), paste0(c("tmax","trange","prcp","cloud","humid","wind"), '.post'))
setkey(c.post, mmsa)

## merge post
post <- merge(r.post, c.post)
rm(r.post, c.post)

## merge
setkey(pre, mmsa, state, lat, lon)
setkey(post, mmsa, state, lat, lon)
diff <- merge(pre, post) # merge together

## calculate total number of years that enter sample
diff[, years := years.pre + years.post]
diff[, subs := subs.pre + subs.post]

## list of variables to difference
v <- gsub('.pre', '', grep('.pre', names(pre), value=TRUE))
invisible(
    foreach::foreach (l=v) %do% {
        diff <- diff[, paste0(l, '.diff') := get(paste0(l, '.post')) - get(paste0(l, '.pre'))]
    } # generate all differences variables
)

diff[, c('years.diff','subs.diff') := NULL]
diff <- diff[!is.nan(tmax.diff),] # remove obs with NaN

ld.sum <- felm(any_ment.diff ~ tmax.diff + prcp.diff + trange.diff + cloud.diff + humid.diff + wind.diff | 0 | 0 | state, data=diff, weights=diff[, subs])
ld.full.sum <- felm(any_ment.diff ~ tmax.diff + prcp.diff + trange.diff + cloud.diff + humid.diff + wind.diff | 0 | 0 | state, data=diff[years==10,], weights=diff[years==10, subs]) # cities with all 10 years in sample
ld.partial.sum <- felm(any_ment.diff ~ tmax.diff + prcp.diff + trange.diff + cloud.diff + humid.diff + wind.diff | 0 | 0 | state, data=diff[years<10,], weights=diff[years<10, subs]) # cities will less than 10 years in sample

# long-diff fall----
pre.years <- 2002:2006
post.years <- 2007:2011

## respondent pre-period averages
r.pre <- brfss[(month(time) %in% c(9,10,11)) & (year(time) %in% pre.years),
               .(any_ment = unlist(lapply(.SD, mean, na.rm=T)),
                 years = uniqueN(year(time)),
                 subs = uniqueN(.I)),
               by=.(mmsa, state, lat, lon),
               .SDcols = c("any_ment")]
setnames(r.pre, c("any_ment", "years", "subs"), paste0(c("any_ment", "years", "subs"), '.pre'))
setkey(r.pre, mmsa)

## city pre-period averages
prism[, trange := tmax - tmin]
c.pre <- prism[(month(time) %in% c(9,10,11)) & (year(time) %in% pre.years),
               lapply(.SD, mean, na.rm=T),
               by=.(mmsa),
               .SDcols = c("tmax","trange","prcp","cloud","humid","wind")]
setnames(c.pre, c("tmax","trange","prcp","cloud","humid","wind"), paste0(c("tmax","trange","prcp","cloud","humid","wind"), '.pre'))
setkey(c.pre, mmsa)

## merge pre
pre <- merge(r.pre, c.pre)
rm(r.pre, c.pre)

## respondent post-period averages
r.post <- brfss[(month(time) %in% c(9,10,11)) & (year(time) %in% post.years),
                .(any_ment = unlist(lapply(.SD, mean, na.rm=T)),
                  years = uniqueN(year(time)),
                  subs = uniqueN(.I)),
                by=.(mmsa, state, lat, lon),
                .SDcols = c("any_ment")]
setnames(r.post, c("any_ment","years","subs"), paste0(c("any_ment","years","subs"), '.post'))
setkey(r.post, mmsa)

## city post-period averages
c.post <- prism[(month(time) %in% c(9,10,11)) & (year(time) %in% post.years),
                lapply(.SD, mean, na.rm=T),
                by=.(mmsa),
                .SDcols = c("tmax","trange","prcp","cloud","humid","wind")]
setnames(c.post, c("tmax","trange","prcp","cloud","humid","wind"), paste0(c("tmax","trange","prcp","cloud","humid","wind"), '.post'))
setkey(c.post, mmsa)

## merge post
post <- merge(r.post, c.post)
rm(r.post, c.post)

## merge
setkey(pre, mmsa, state, lat, lon)
setkey(post, mmsa, state, lat, lon)
diff <- merge(pre, post) # merge together

## calculate total number of years that enter sample
diff[, years := years.pre + years.post]
diff[, subs := subs.pre + subs.post]

## list of variables to difference
v <- gsub('.pre', '', grep('.pre', names(pre), value=TRUE))
invisible(
    foreach::foreach (l=v) %do% {
        diff <- diff[, paste0(l, '.diff') := get(paste0(l, '.post')) - get(paste0(l, '.pre'))]
    } # generate all differences variables
)

diff[, c('years.diff','subs.diff') := NULL]
diff <- diff[!is.nan(tmax.diff),] # remove obs with NaN

ld.fall <- felm(any_ment.diff ~ tmax.diff + prcp.diff + trange.diff + cloud.diff + humid.diff + wind.diff | 0 | 0 | state, data=diff, weights=diff[, subs])
ld.full.fall <- felm(any_ment.diff ~ tmax.diff + prcp.diff + trange.diff + cloud.diff + humid.diff + wind.diff | 0 | 0 | state, data=diff[years==10,], weights=diff[years==10, subs]) # cities with all 10 years in sample
ld.partial.fall <- felm(any_ment.diff ~ tmax.diff + prcp.diff + trange.diff + cloud.diff + humid.diff + wind.diff | 0 | 0 | state, data=diff[years<10,], weights=diff[years<10, subs]) # cities will less than 10 years in sample

# long-diff winter----
pre.years <- 2002:2006
post.years <- 2007:2011

## respondent pre-period averages
r.pre <- brfss[(month(time) %in% c(12,1,2)) & (year(time) %in% pre.years),
               .(any_ment = unlist(lapply(.SD, mean, na.rm=T)),
                 years = uniqueN(year(time)),
                 subs = uniqueN(.I)),
               by=.(mmsa, state, lat, lon),
               .SDcols = c("any_ment")]
setnames(r.pre, c("any_ment", "years", "subs"), paste0(c("any_ment", "years", "subs"), '.pre'))
setkey(r.pre, mmsa)

## city pre-period averages
prism[, trange := tmax - tmin]
c.pre <- prism[(month(time) %in% c(12,1,2)) & (year(time) %in% pre.years),
               lapply(.SD, mean, na.rm=T),
               by=.(mmsa),
               .SDcols = c("tmax","trange","prcp","cloud","humid","wind")]
setnames(c.pre, c("tmax","trange","prcp","cloud","humid","wind"), paste0(c("tmax","trange","prcp","cloud","humid","wind"), '.pre'))
setkey(c.pre, mmsa)

## merge pre
pre <- merge(r.pre, c.pre)
rm(r.pre, c.pre)

## respondent post-period averages
r.post <- brfss[(month(time) %in% c(12,1,2)) & (year(time) %in% post.years),
                .(any_ment = unlist(lapply(.SD, mean, na.rm=T)),
                  years = uniqueN(year(time)),
                  subs = uniqueN(.I)),
                by=.(mmsa, state, lat, lon),
                .SDcols = c("any_ment")]
setnames(r.post, c("any_ment","years","subs"), paste0(c("any_ment","years","subs"), '.post'))
setkey(r.post, mmsa)

## city post-period averages
c.post <- prism[(month(time) %in% c(12,1,2)) & (year(time) %in% post.years),
                lapply(.SD, mean, na.rm=T),
                by=.(mmsa),
                .SDcols = c("tmax","trange","prcp","cloud","humid","wind")]
setnames(c.post, c("tmax","trange","prcp","cloud","humid","wind"), paste0(c("tmax","trange","prcp","cloud","humid","wind"), '.post'))
setkey(c.post, mmsa)

## merge post
post <- merge(r.post, c.post)
rm(r.post, c.post)

## merge
setkey(pre, mmsa, state, lat, lon)
setkey(post, mmsa, state, lat, lon)
diff <- merge(pre, post) # merge together

## calculate total number of years that enter sample
diff[, years := years.pre + years.post]
diff[, subs := subs.pre + subs.post]

## list of variables to difference
v <- gsub('.pre', '', grep('.pre', names(pre), value=TRUE))
invisible(
    foreach::foreach (l=v) %do% {
        diff <- diff[, paste0(l, '.diff') := get(paste0(l, '.post')) - get(paste0(l, '.pre'))]
    } # generate all differences variables
)

diff[, c('years.diff','subs.diff') := NULL]
diff <- diff[!is.nan(tmax.diff),] # remove obs with NaN

ld.win <- felm(any_ment.diff ~ tmax.diff + prcp.diff + trange.diff + cloud.diff + humid.diff + wind.diff | 0 | 0 | state, data=diff, weights=diff[, subs])
ld.full.win <- felm(any_ment.diff ~ tmax.diff + prcp.diff + trange.diff + cloud.diff + humid.diff + wind.diff | 0 | 0 | state, data=diff[years==10,], weights=diff[years==10, subs]) # cities with all 10 years in sample
ld.partial.win <- felm(any_ment.diff ~ tmax.diff + prcp.diff + trange.diff + cloud.diff + humid.diff + wind.diff | 0 | 0 | state, data=diff[years<10,], weights=diff[years<10, subs]) # cities will less than 10 years in sample
# katrina diff-n-diff regression----

## set up treatment variables
brfss[, treat := (state %like% "LA|MS") | (name %like% "Birmingham AL|Mobile AL|Tuscaloosa AL|Miami FL|Key West FL|Naples FL|Panama City FL|Pensacola FL")] 

## set pre-Katrina as pre and post-Katrina as post
brfss[, post := time > "2005-08-24"]

## interaction
brfss[, postXtreat := post*treat]

katrinareg <- summary(felm(any_ment ~ post + treat + postXtreat
             | 0 
             | 0 
             | state, 
             data=brfss, 
             exactDOF=TRUE))
```

# Figure 1

```{r fig1, eval=FALSE}
# future hot days----
## trim projection rasters to remove outer NA values
year2010 <- trim(bcsd.tmax[[1]], padding=0, values=NA)
year2050 <- trim(bcsd.tmax[[2]], padding=0, values=NA)
year2099 <- trim(bcsd.tmax[[3]], padding=0, values=NA)

## stack these together
forecast <- stack(year2010, year2050, year2099)
id <- c('2010','2050','2099') # set year id
forecast <- setZ(forecast, id) # set year as z variable
names(forecast) <- c('2010','2050','2099') # set names

## plot raster forecasts
theme = rasterTheme(region=c('gray95','gray60','#fc8d62','red4')) # define color palette

p <- levelplot(forecast, 
               par.settings = theme, #set theme
               xlab=NULL, # kill x-axis
               ylab=NULL, # kill y-axis
               col=NA,
               margin=TRUE, # remove margin padding
               scales=list(draw=FALSE),
               layout=c(3,1), # horizontal layout
               names.attr=c("", "", ""), # inset titles
               colorkey=list(space="right", labels=list(cex=2.5), height=0.75), # make legend on bottom of plot
               panel=panel.levelplot.raster, # makes it so that grid lines are not plotted
               pretty=TRUE,
               par.strip.text=list(cex=2.5, lines=1.5),
               sub=list("Days with Maxmimum Temperature >30°C", cex=3)
               
    )
p <- update(p, par.settings = 
                list(axis.line = list(col = 0), 
                     panel.background=list(col="transparent")
                     )
            ) # remove plot borders and add subtitle below legend
p <- p + layer(sp.lines(states, lwd=0.5, col="gray60")) # add us shapefiles

# future prcp days----
## trim forecast rasters to remove outer NA values
year2010 <- trim(bcsd.prcp.days[[1]], padding=0, values=NA)
year2050 <- trim(bcsd.prcp.days[[2]], padding=0, values=NA)
year2099 <- trim(bcsd.prcp.days[[3]], padding=0, values=NA)

## stack these together
forecast <- stack(year2010, year2050, year2099)
id <- c('2010','2050','2099') # set year id
forecast <- setZ(forecast, id) # set year as z variable
names(forecast) <- c('2010','2050','2099') # set names

## plot raster forecasts
theme = rasterTheme(region=c('gray95','gray60','#8da0cb','blue4')) # define color palette

p1 <- levelplot(forecast, 
               par.settings = theme, #set theme
               xlab=NULL, # kill x-axis
               ylab=NULL, # kill y-axis
               col=NA,
               margin=TRUE, # remove margin padding
               scales=list(draw=FALSE),
               layout=c(3,1), # horizontal layout
               names.attr=c("", "", ""), # inset titles
               colorkey=list(space="right", labels=list(cex=2.5), height=0.75), # make legend on bottom of plot
               panel=panel.levelplot.raster, # makes it so that grid lines are not plotted
               pretty=TRUE,
               par.strip.text=list(cex=2.5, lines=1.5),
               sub=list("Days with Measurable Precipitation", cex=3)
               
    )
p1 <- update(p1, par.settings = 
                list(axis.line = list(col = 0), 
                     panel.background=list(col="transparent")
                     )
            ) # remove plot borders and add subtitle below legend
p1 <- p1 + layer(sp.lines(states, lwd=0.5, col="gray60")) # add us shapefiles

# plot of marginal effects from main model----
cut <- felm(any_ment ~ cutprcp + cuttmax + avg.trange + avg.humid + avg.cloud + avg.wind
          | time + mmsa + stateyr
          | 0 
          | state
          , data=brfss,
          exactDOF=TRUE)
cut.plot <- bin.plot(felm.est = cut,
                     plotvar = "tmax",
                     breaks = 5,
                     omit = c(10,15),
                     ylimit = c(-0.75, 3),
                     xlabel = "30-Day Avg. Max. Temp. in \u2103",
                     ylabel = "P(Report Mental Health Issue)\n(Percentage Points)",
                     y.axis.title = element_blank(),
                     y.axis.text = element_blank(),
                     y.axis.line = element_blank(),
                     y.axis.ticks = element_blank(),
                     panel = ""
                     )
hist <- axis_canvas(cut.plot, axis = "x") + 
        geom_histogram(data = prism, aes(x = avg.tmax), bins=30, fill="#8da0cb", color='gray60', alpha=0.5, size=0.1)
cut.plot <- insert_xaxis_grob(cut.plot, hist, position = "bottom", height = grid::unit(0.05, "null"))
cut.plot <- ggdraw(cut.plot)

# prcp days----
prcp <- bin.plot(felm.est = cut,
                     plotvar = "prcp",
                     breaks = 5,
                     omit = c(0,0),
                     ylimit = c(-0.75, 3),
                     xlabel = "Days of Precipitation",
                     ylabel = "P(Report Mental Health Issue)\n(Percentage Points)",
                     y.axis.title = element_blank(),
                     y.axis.text = element_blank(),
                     y.axis.line = element_blank(),
                     y.axis.ticks = element_blank(),
                     panel = ""
                     )
hist <- axis_canvas(prcp, axis = "x") + 
        geom_histogram(data = prism, aes(x = prcp.days), bins=31, fill="#8da0cb", color='gray60', alpha=0.5, size=0.1)
prcp <- insert_xaxis_grob(prcp, hist, position = "bottom", height = grid::unit(0.05, "null"))
prcp <- ggdraw(prcp)

# just y-label----
ylabel <- bin.plot(felm.est = cut,
                     plotvar = "tmax",
                     breaks = 5,
                     omit = c(10,15),
                     ylimit = c(-0.75, 3),
                     xlabel = "30-Day Avg. Max. Temp. in \u2103",
                     ylabel = "P(Mental Health Issue)\n(Percentage Points)",
                     panel = ""
                     )
hist <- axis_canvas(ylabel, axis = "x") + 
        geom_histogram(data = prism, aes(x = tmax), bins=30, fill="#8da0cb", color='gray60', alpha=0.5, size=0.1)
ylabel <- insert_xaxis_grob(ylabel, hist, position = "bottom", height = grid::unit(0.05, "null"))
ylabel <- gtable_filter(ylabel, 'axis-l|ylab', trim=F)
ylabel <- ggdraw(ylabel)

# plot----
cairo_pdf(filename = './figure1.pdf', width=13.5, height=12.5, fallback_resolution = 1200)
grid.arrange( 
             p,
             p1,
             arrangeGrob(rectGrob(gp=gpar(col="white")), arrangeGrob(ylabel, cut.plot, prcp, ncol=3, widths=c(0.08, 0.46, 0.46)), rectGrob(gp=gpar(col="white")), ncol=3, widths=c(0.04, 0.92, 0.04)),
             nrow=3, ncol=1, heights=c(0.3, 0.3, 0.4))
grid.text("a", x=unit(0.02, "npc"), y=unit(0.96, "npc"), gp=gpar(fontsize=40,font=2))
grid.text("b", x=unit(0.02, "npc"), y=unit(0.65, "npc"), gp=gpar(fontsize=40,font=2))
grid.text("c", x=unit(0.14, "npc"), y=unit(0.1, "npc"), gp=gpar(fontsize=40,font=2))
grid.text("d", x=unit(0.55, "npc"), y=unit(0.1, "npc"), gp=gpar(fontsize=40,font=2))
grid.text("2010", x=unit(0.15, "npc"), y=unit(0.85, "npc"), gp=gpar(fontsize=40,font=2, col="white"))
grid.text("2050", x=unit(0.45, "npc"), y=unit(0.85, "npc"), gp=gpar(fontsize=40,font=2, col="white"))
grid.text("2099", x=unit(0.75, "npc"), y=unit(0.85, "npc"), gp=gpar(fontsize=40,font=2, col="white"))
grid.text("2010", x=unit(0.15, "npc"), y=unit(0.55, "npc"), gp=gpar(fontsize=40,font=2, col="white"))
grid.text("2050", x=unit(0.45, "npc"), y=unit(0.55, "npc"), gp=gpar(fontsize=40,font=2, col="white"))
grid.text("2099", x=unit(0.75, "npc"), y=unit(0.55, "npc"), gp=gpar(fontsize=40,font=2, col="white"))
dev.off()
```

```{r fig.contemporaneous, out.width="100%", fig.cap="Higher temperatures and precipitation increase risk of mental health issues. Panel (a) presents projections of warming in the US along the RCP8.5 high emissions scenario. The annual number of days with maximum temperatures >30$^\\circ{}$C is projected to markedly increase over this century. Panel (b) plots the projected number of annual days with measurable precipitation (>1mm). Panel (c) draws from nearly 2 million respondents' reports of monthly mental health issues between 2002 and 2012. It plots the predicted probability of reporting any mental health issues for each thirty-day average maximum temperature bin. The probability of mental health issues steadily increases past 10$^\\circ{}$C-15$^\\circ{}$C. Panel (d) depicts that precipitation days increase the probability of mental health issues. Dotted lines represent 95\\% confidence intervals."}
knitr::include_graphics("figure1.pdf")
```

```{r}
contemp <- fig_num(name = "contemp", caption = "")
```

# Figure 2

```{r fig2, eval=FALSE}
# first versus fourth quartile income----
poor <- felm(any_ment ~ cutprcp + cuttmax + avg.trange + avg.humid + avg.cloud + avg.wind
            | time + mmsa + stateyr
            | 0
            | state
            , data=brfss[income2 <= 22500,],
            exactDOF=TRUE)

rich <- felm(any_ment ~ cutprcp + cuttmax + avg.trange + avg.humid + avg.cloud + avg.wind
            | time + mmsa + stateyr
            | 0
            | state
            , data=brfss[income2 >= 90000,],
            exactDOF=TRUE)

# get coefficients of variables
beta.poor <- poor$coefficients[row.names(poor$coefficients)=="cuttmax(30, Inf]"]
beta.rich <- rich$coefficients[row.names(rich$coefficients)=="cuttmax(30, Inf]"]

# cluster robust standard errors
se.poor <- poor$cse[names(poor$cse)=="cuttmax(30, Inf]"]
se.rich <- rich$cse[names(rich$cse)=="cuttmax(30, Inf]"]

# upper and lower confidence bounds
upper.poor <- beta.poor + se.poor
lower.poor <- beta.poor - se.poor
upper.rich <- beta.rich + se.rich
lower.rich <- beta.rich - se.rich

d <- as.data.table(cbind(rbind(beta.poor,beta.rich),rbind(upper.poor,upper.rich),rbind(lower.poor,lower.rich),rbind("Low","High")))
setnames(d, c('beta', 'upper_bound', 'lower_bound', 'variables'))
d <- d[, beta := as.numeric(beta)]
d <- d[, upper_bound := as.numeric(upper_bound)]
d <- d[, lower_bound := as.numeric(lower_bound)]
d <- d[, variables := factor(variables, levels=c("Low", "High"))]

inc <- ggplot(data = d, aes(variables, beta, fill=variables)) +
        geom_bar(stat = "identity", aes(width = 0.8), color='black') +
        geom_linerange(aes(ymax = upper_bound, ymin = lower_bound)) +
        scale_fill_manual(values=c("#fc8d62","gray70"),breaks=c("Low","High")) +
        coord_cartesian(xlim = c(0.5, 2.5), ylim = c(0, 3.1), expand=FALSE) +
        ylab("Marginal Effect of >30\u2103 on\nP(Mental Health Issue)\n(Percentage Points)") +
        xlab("Income") +
        geom_hline(yintercept = 0) +
        theme(title = element_text(size=20, angle=0, vjust=.5, face="bold"),
              plot.title = element_text(hjust=0.025),
              axis.title.x = element_text(size=40, angle=0, vjust=-0, face="bold"),
              axis.title.y = element_blank(),
              axis.text.y = element_blank(),
              axis.text.x = element_text(size=40, face="bold", color="gray30"),
              panel.grid.major = element_blank(),
              panel.grid.minor = element_blank(),
              panel.background = element_blank(),
              axis.line.x = element_line(),
              axis.line.y = element_blank(),
              axis.ticks.y = element_blank(),
              legend.position="none"
        )

# female vs male----
female <- felm(any_ment ~ cutprcp + cuttmax + avg.trange + avg.humid + avg.cloud + avg.wind
            | time + mmsa + stateyr
            | 0
            | state
            , data=brfss[female == 1,],
            exactDOF=TRUE)

male <- felm(any_ment ~ cutprcp + cuttmax + avg.trange + avg.humid + avg.cloud + avg.wind
            | time + mmsa + stateyr
            | 0
            | state
            , data=brfss[female == 0,],
            exactDOF=TRUE)

# get coefficients of variables
beta.female <- female$coefficients[row.names(female$coefficients)=="cuttmax(30, Inf]"]
beta.male <- male$coefficients[row.names(male$coefficients)=="cuttmax(30, Inf]"]

# cluster robust standard errors
se.female <- female$cse[names(female$cse)=="cuttmax(30, Inf]"]
se.male <- male$cse[names(male$cse)=="cuttmax(30, Inf]"]

# upper and lower confidence bounds
upper.female <- beta.female + se.female
lower.female <- beta.female - se.female
upper.male <- beta.male + se.male
lower.male <- beta.male - se.male

d <- as.data.table(cbind(rbind(beta.female,beta.male),rbind(upper.female,upper.male),rbind(lower.female,lower.male),rbind("Women","Men")))
setnames(d, c('beta', 'upper_bound', 'lower_bound', 'variables'))
d <- d[, beta := as.numeric(beta)]
d <- d[, upper_bound := as.numeric(upper_bound)]
d <- d[, lower_bound := as.numeric(lower_bound)]
d <- d[, variables := factor(variables, levels=c("Women", "Men"))]

gender <- ggplot(data = d, aes(variables, beta, fill=variables)) +
        geom_bar(stat = "identity", aes(width = 0.8), color='black') +
        geom_linerange(aes(ymax = upper_bound, ymin = lower_bound)) +
        scale_fill_manual(values=c("#fc8d62","gray70"),breaks=c("Women","Men")) +
        coord_cartesian(xlim = c(0.5, 2.5), ylim = c(0, 3.1), expand=FALSE) +
        ylab("Marginal Effect of >30\u2103 on\nP(Mental Health Issue)\n(Percentage Points)") +
        xlab("Gender") +
        geom_hline(yintercept = 0) +
        theme(title = element_text(size=20, angle=0, vjust=.5, face="bold"),
              plot.title = element_text(hjust=0.025),
              axis.title.x = element_text(size=40, angle=0, vjust=-0, face="bold"),
              axis.title.y = element_blank(),
              axis.text.y = element_blank(),
              axis.text.x = element_text(size=40, face="bold", color="gray30"),
              panel.grid.major = element_blank(),
              panel.grid.minor = element_blank(),
              panel.background = element_blank(),
              axis.line.x = element_line(),
              axis.line.y = element_blank(),
              axis.ticks.y = element_blank(),
              legend.position="none"
        )

# justy ----
justy <- ggplot(data = d, aes(variables, beta, fill=variables)) +
        geom_bar(stat = "identity", aes(width = 0.8), color='black') +
        geom_linerange(aes(ymax = upper_bound, ymin = lower_bound)) +
        scale_fill_manual(values=c("#fc8d62","gray70"),breaks=c("Female","Male")) +
        coord_cartesian(xlim = c(0.5, 2.5), ylim = c(0, 3.1), expand=FALSE) +
        ylab("Marginal Effect of >30\u2103 on\nP(Mental Health Issue)\n(Percentage Points)") +
        xlab("Gender") +
        geom_hline(yintercept = 0) +
        theme(title = element_text(size=20, angle=0, vjust=.5, face="bold"),
              plot.title = element_text(hjust=0.025),
              axis.title.x = element_text(size=40, angle=0, vjust=-0, face="bold"),
              axis.title.y = element_text(size=30, angle=90, vjust=1, face="bold"),
              axis.text.y = element_text(size=30, angle=0, vjust=.5, face="bold"),
              axis.text.x = element_text(size=40, face="bold", color="gray30"),
              panel.grid.major = element_blank(),
              panel.grid.minor = element_blank(),
              panel.background = element_blank(),
              axis.line.y = element_line(),
              axis.line.x = element_line(),
              legend.position="none"
        )
justy <- gtable_filter(ggplotGrob(justy), 'axis-l|ylab', trim=F)
justy <- ggdraw(justy)

# plot----
cairo_pdf(filename = "./figure2.pdf", width=17, height=6.5, fallback_resolution = 1200)
grid.arrange(justy, inc, gender, ncol=3, widths=c(0.1, 0.45, 0.45))
grid.text("a", x=unit(0.12, "npc"), y=unit(0.97, "npc"), gp=gpar(fontsize=40,font=2))
grid.text("b", x=unit(0.58, "npc"), y=unit(0.97, "npc"), gp=gpar(fontsize=40,font=2))
dev.off()
```

```{r fig.hetero, out.width="100%", fig.cap="Heat effects are more acute among low income respondents and among women. Panel (a) displays the marginal effects produced by splitting the sample between the first and fourth quartiles of income. Those with incomes in the first quartile have amplified probabilities of reporting mental health issues in response to higher temperatures. Panel (b) stratifies the sample by gender. The effects of higher temperatures on mental health issues are larger among women. All marginal effects are significant at the $\\alpha = 0.05$ level. Error bars are SEM."}
knitr::include_graphics("figure2.pdf")
```

```{r}
hetero <- fig_num(name = "hetero", caption = "")
```

# Figure 3

```{r fig3, eval=FALSE}
# long diff----
pre.years <- 2002:2006
post.years <- 2007:2011

## respondent pre-period averages
r.pre <- brfss[year(time) %in% pre.years,
               .(any_ment = unlist(lapply(.SD, mean, na.rm=T)),
                 years = uniqueN(year(time)),
                 subs = uniqueN(.I)),
               by=.(mmsa, state, lat, lon),
               .SDcols = c("any_ment")]
setnames(r.pre, c("any_ment", "years", "subs"), paste0(c("any_ment", "years", "subs"), '.pre'))
setkey(r.pre, mmsa)

## city pre-period averages
prism[, trange := tmax - tmin]
c.pre <- prism[year(time) %in% pre.years,
               lapply(.SD, mean, na.rm=T),
               by=.(mmsa),
               .SDcols = c("tmax","trange","prcp","cloud","humid","wind")]
setnames(c.pre, c("tmax","trange","prcp","cloud","humid","wind"), paste0(c("tmax","trange","prcp","cloud","humid","wind"), '.pre'))
setkey(c.pre, mmsa)

## merge pre
pre <- merge(r.pre, c.pre)
rm(r.pre, c.pre)

## respondent post-period averages
r.post <- brfss[year(time) %in% post.years,
                .(any_ment = unlist(lapply(.SD, mean, na.rm=T)),
                  years = uniqueN(year(time)),
                  subs = uniqueN(.I)),
                by=.(mmsa, state, lat, lon),
                .SDcols = c("any_ment")]
setnames(r.post, c("any_ment","years","subs"), paste0(c("any_ment","years","subs"), '.post'))
setkey(r.post, mmsa)

## city post-period averages
c.post <- prism[year(time) %in% post.years,
                lapply(.SD, mean, na.rm=T),
                by=.(mmsa),
                .SDcols = c("tmax","trange","prcp","cloud","humid","wind")]
setnames(c.post, c("tmax","trange","prcp","cloud","humid","wind"), paste0(c("tmax","trange","prcp","cloud","humid","wind"), '.post'))
setkey(c.post, mmsa)

## merge post
post <- merge(r.post, c.post)
rm(r.post, c.post)

## merge
setkey(pre, mmsa, state, lat, lon)
setkey(post, mmsa, state, lat, lon)
diff <- merge(pre, post) # merge together

## calculate total number of years that enter sample
diff[, years := years.pre + years.post]
diff[, subs := subs.pre + subs.post]

## list of variables to difference
v <- gsub('.pre', '', grep('.pre', names(pre), value=TRUE))
invisible(
    foreach::foreach (l=v) %do% {
        diff <- diff[, paste0(l, '.diff') := get(paste0(l, '.post')) - get(paste0(l, '.pre'))]
    } # generate all differences variables
)

diff[, c('years.diff','subs.diff') := NULL]
diff <- diff[!is.nan(tmax.diff),] # remove obs with NaN

ld <- felm(any_ment.diff ~ tmax.diff + prcp.diff + trange.diff + cloud.diff + humid.diff + wind.diff | 0 | 0 | state, data=diff, weights=diff[, subs])
ld.full <- felm(any_ment.diff ~ tmax.diff + prcp.diff + trange.diff + cloud.diff + humid.diff + wind.diff | 0 | 0 | state, data=diff[years==10,], weights=diff[years==10, subs]) # cities with all 10 years in sample
ld.partial <- felm(any_ment.diff ~ tmax.diff + prcp.diff + trange.diff + cloud.diff + humid.diff + wind.diff | 0 | 0 | state, data=diff[years<10,], weights=diff[years<10, subs]) # cities will less than 10 years in sample

# calculate all possible combinations of long-difference years----
# a <- foreach (i=1:10, .combine = rbind) %dopar% {
#         foreach (j=split(combn(2002:2011, i), col(combn(2002:2011, i))), .combine = rbind) %do% {
#             foreach (k=1:10, .combine = rbind) %do% {
#                 foreach (m=split(combn(2002:2011, k), col(combn(2002:2011, k))), .combine = rbind) %do% {
#                         pre.years <- j
#                         post.years <- m
#                         if (all(post.years - pre.years >= 5) & all(!(pre.years %in% post.years))) { # ensure that years do not overlap and that post is strictly greatereq 5 years from pre-period
# 
#                             ## respondent pre-period averages
#                             r.pre <- brfss[year(time) %in% pre.years,
#                                            .(any_ment = unlist(lapply(.SD, mean, na.rm=T)),
#                                              years = uniqueN(year(time)),
#                                              subs = uniqueN(.I)),
#                                            by=.(mmsa, state, lat, lon),
#                                            .SDcols = c("any_ment")]
#                             setnames(r.pre, c("any_ment", "years", "subs"), paste0(c("any_ment", "years", "subs"), '.pre'))
#                             setkey(r.pre, mmsa)
# 
#                             ## city pre-period averages
#                             prism[, trange := tmax - tmin]
#                             c.pre <- prism[year(time) %in% pre.years,
#                                            lapply(.SD, mean, na.rm=T),
#                                            by=.(mmsa),
#                                            .SDcols = c("tmax","trange","prcp","cloud","humid","wind")]
#                             setnames(c.pre, c("tmax","trange","prcp","cloud","humid","wind"), paste0(c("tmax","trange","prcp","cloud","humid","wind"), '.pre'))
#                             setkey(c.pre, mmsa)
# 
#                             ## merge pre
#                             pre <- merge(r.pre, c.pre)
#                             rm(r.pre, c.pre)
# 
#                             ## respondent post-period averages
#                             r.post <- brfss[year(time) %in% post.years,
#                                             .(any_ment = unlist(lapply(.SD, mean, na.rm=T)),
#                                               years = uniqueN(year(time)),
#                                               subs = uniqueN(.I)),
#                                             by=.(mmsa, state, lat, lon),
#                                             .SDcols = c("any_ment")]
#                             setnames(r.post, c("any_ment","years","subs"), paste0(c("any_ment","years","subs"), '.post'))
#                             setkey(r.post, mmsa)
# 
#                             ## city post-period averages
#                             c.post <- prism[year(time) %in% post.years,
#                                             lapply(.SD, mean, na.rm=T),
#                                             by=.(mmsa),
#                                             .SDcols = c("tmax","trange","prcp","cloud","humid","wind")]
#                             setnames(c.post, c("tmax","trange","prcp","cloud","humid","wind"), paste0(c("tmax","trange","prcp","cloud","humid","wind"), '.post'))
#                             setkey(c.post, mmsa)
# 
#                             ## merge post
#                             post <- merge(r.post, c.post)
#                             rm(r.post, c.post)
# 
#                             ## merge
#                             setkey(pre, mmsa, state, lat, lon)
#                             setkey(post, mmsa, state, lat, lon)
#                             diff <- merge(pre, post) # merge together
# 
#                             ## calculate total number of years that enter sample
#                             diff[, years := years.pre + years.post]
#                             diff[, subs := subs.pre + subs.post]
# 
#                             ## list of variables to difference
#                             v <- gsub('.pre', '', grep('.pre', names(pre), value=TRUE))
#                             invisible(
#                                 foreach::foreach (l=v) %do% {
#                                     diff <- diff[, paste0(l, '.diff') := get(paste0(l, '.post')) - get(paste0(l, '.pre'))]
#                                 } # generate all differences variables
#                             )
#                             rm(pre, post, l, v)
# 
#                             diff[, c('years.diff','subs.diff') := NULL]
#                             diff <- diff[!is.nan(tmax.diff),] # remove obs with NaN
# 
#                             ld.1 <- felm(any_ment.diff ~ tmax.diff + prcp.diff + trange.diff + cloud.diff + humid.diff + wind.diff | 0 | 0 | state, data=diff, weights=diff[, subs])
# 
#                             d <- data.table(coef = ld.1$coef[2], cse = ld.1$cse[2], post = paste(m, collapse=""), pre = paste(j, collapse=""))
#                         } else {
#                             d <- data.table(coef = NA, cse = NA, post = paste(m, collapse=""), pre = paste(j, collapse=""))
#                         }
#                         print(j)
#                         print(m)
#                         return(d)
#                 }
#             }
#         }
# }
# 
# saveRDS(a, "./long_diff_combinations.rds", compress=TRUE)

# long diff all possible combinations of post-pre----

## calculate significant parameters
ldcomb[, upr := coef + 1.96*cse]
ldcomb[, lwr := coef - 1.96*cse]
ldcomb[, sig := upr < 0 | lwr > 0]
ldcomb[, hist.coef := coef * -1] # ggplot hack to get this to display properly

## select significant parameters and order by magnitude
ld.sig <- ldcomb[sig==TRUE,]
setkey(ld.sig, coef)
ld.sig[, iter := .I]

# plots----
## difference by location
us <- map_data("state")
map <- ggplot(data=us, aes(x=long, y=lat, group=group)) + 
            geom_path(color='gray40') +
            geom_point(data=diff, aes(x=lon, y=lat, group=NULL, size=subs, fill=any_ment.diff), pch=21, color='white', alpha=0.4) + 
            scale_size_continuous(range = c(5,25)) +
            scale_fill_gradientn(colours=c("blue4","#8da0cb","gray80","#fc8d62","red4")) +
            guides(fill = guide_colorbar(barwidth = 23, barheight = 1, ticks=FALSE, title.position="top", title.hjust = 0.5, direction = "horizontal"), size = FALSE) +
            labs(fill="\u0394 Mental Health Issues\n(from 2002-2006 to 2007-2011, pp)") +
          theme(
          axis.title.x = element_blank(),
          axis.title.y = element_blank(),
          axis.text.y = element_blank(),
          axis.text.x = element_blank(),
          axis.ticks = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.background = element_blank(),
          axis.line = element_blank(),
          legend.title = element_text(size=25, angle=0, face="bold"),
          legend.position = "none",
          legend.text = element_text(size=20, face="bold")
          )

legend1 <- ggplot(data=us, aes(x=long, y=lat, group=group)) + 
            geom_path(color='gray40') +
            geom_point(data=diff, aes(x=lon, y=lat, group=NULL, size=subs, fill=any_ment.diff), pch=21, color='white', alpha=0.4) + 
            scale_size_continuous(range = c(5,25)) +
            scale_fill_gradientn(colours=c("blue4","#8da0cb","gray80","#fc8d62","red4")) +
            guides(fill = guide_colorbar(barwidth = 23, barheight = 1, ticks=FALSE, title.position="top", title.hjust = 0.5, direction = "horizontal"), size = FALSE) +
            labs(fill="\u0394 Mental Health Issues\n(2002-2006 to 2007-2011, pp)") +
          theme(
          axis.title.x = element_blank(),
          axis.title.y = element_blank(),
          axis.text.y = element_blank(),
          axis.text.x = element_blank(),
          axis.ticks = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.background = element_blank(),
          axis.line = element_blank(),
          legend.title = element_text(size=25, angle=0, face="bold"),
          legend.position = c(0.155,0.7),
          legend.text = element_text(size=20, face="bold")
          )
legend1 <- gtable_filter(ggplotGrob(legend1), 'guide-box', trim=T)
legend1 <- ggdraw(legend1)

## prism differences
pdiff <- as.data.table(as(prism.diff, "SpatialPixelsDataFrame"))
p <- ggplot(data=us, aes(x=long, y=lat, group=group)) + 
            geom_path(color='gray40') +
            geom_tile(data=pdiff, aes(x=x, y=y, group=NULL, fill=layer)) +          
            scale_fill_gradientn(colours=c("blue4","#8da0cb","gray90","#fc8d62","orangered3","red3","orangered4","red4")) +
            guides(fill = guide_colorbar(barwidth = 23, barheight = 1, ticks=FALSE, title.position="top", title.hjust = 0.5, direction = "horizontal"), size = FALSE) +
            labs(fill="\u0394 Avg. Max. Temp.\n(2002-2006 to 2007-2011, \u2103)") +
          theme(
          axis.title.x = element_blank(),
          axis.title.y = element_blank(),
          axis.text.y = element_blank(),
          axis.text.x = element_blank(),
          axis.ticks = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.background = element_blank(),
          axis.line = element_blank(),
          legend.title = element_text(size=25, angle=0, face="bold"),
          legend.position = "none",
          legend.text = element_text(size=20, face="bold")
          )

legend2 <- ggplot(data=us, aes(x=long, y=lat, group=group)) + 
            geom_path(color='gray40') +
            geom_tile(data=pdiff, aes(x=x, y=y, group=NULL, fill=layer)) +          
            scale_fill_gradientn(colours=c("blue4","#8da0cb","gray90","#fc8d62","orangered3","red3","orangered4","red4")) +
            guides(fill = guide_colorbar(barwidth = 23, barheight = 1, ticks=FALSE, title.position="top", title.hjust = 0.5, direction = "horizontal"), size = FALSE) +
            labs(fill="\u0394 Avg. Max. Temp.\n(2002-2006 to 2007-2011, \u2103)") +
          theme(
          axis.title.x = element_blank(),
          axis.title.y = element_blank(),
          axis.text.y = element_blank(),
          axis.text.x = element_blank(),
          axis.ticks = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.background = element_blank(),
          axis.line = element_blank(),
          legend.title = element_text(size=25, angle=0, face="bold"),
          legend.position = c(0.155,0.7),
          legend.text = element_text(size=20, face="bold")
          )
legend2 <- gtable_filter(ggplotGrob(legend2), 'guide-box', trim=T)
legend2 <- legend <- ggdraw(legend2)

## preferred specification estimates
ld <- data.table(coef = c(ld$coef[2], ld.full$coef[2], ld.partial$coef[2]), se = c(ld$cse[2], ld.full$cse[2], ld.partial$cse[2]), model = c("All", "Full", "Missing"))
ld[, upr := coef + 1.96*se]
ld[, lwr := coef - 1.96*se]
coef <- ggplot(ld, aes(model, coef)) +
        geom_hline(yintercept=0, lty=2, lwd=1, colour="grey60") +
        geom_errorbar(aes(ymin=lwr, ymax=upr), lwd=2, colour="#8da0cb", width=0) +
        geom_errorbar(aes(ymin=coef - se, ymax=coef + se), lwd=3.5, colour="#fc8d62", width=0) +
        geom_point(size=4, pch=21, fill="red4") +
        ylab("+1\u0394\u2103 on \u0394 Mental Health Issues\n(Percentage Points)") +
        xlab("Long-Diff. Model") +
        coord_cartesian(ylim=c(-2.5,5)) +
        theme(plot.title = element_text(face="bold", size=40),
              axis.title.x = element_text(size=30, angle=0, face="bold"),
              axis.title.y = element_blank(),
              axis.text.y = element_blank(),
              axis.text.x = element_text(size=25, face="bold"),
              panel.grid.major = element_blank(),
              panel.grid.minor = element_blank(),
              panel.background = element_blank(),
              axis.line.x = element_line(),
              axis.line.y = element_blank(),
              axis.ticks.y = element_blank(),
              legend.title = element_blank()
        )

## justy
justy <- ggplot(ld, aes(model, coef)) +
        geom_hline(yintercept=0, lty=2, lwd=1, colour="grey60") +
        geom_errorbar(aes(ymin=lwr, ymax=upr), lwd=2, colour="#8da0cb", width=0) +
        geom_errorbar(aes(ymin=coef - se, ymax=coef + se), lwd=5, colour="#fc8d62", width=0) +
        geom_point(size=8, pch=21, fill="red4") +
        ylab("Effect of +1\u0394\u2103 on \n\u0394 Mental Health Issues\n(Percentage Points)") +
        xlab("Long-Diff. Model") +
        coord_cartesian(ylim=c(-2.5,5)) +
        theme(plot.title = element_text(face="bold", size=40),
              axis.title.x = element_text(size=30, angle=0, face="bold"),
              axis.title.y = element_text(size=25, angle=90, vjust=-0.25, face="bold"),
              axis.text.y = element_text(size=20, angle=0, face="bold") ,
              axis.text.x = element_text(size=25, face="bold"),
              panel.grid.major = element_blank(),
              panel.grid.minor = element_blank(),
              panel.background = element_blank(),
              axis.line.x = element_line(),
              axis.line.y = element_line(),
              legend.title = element_blank()
        )
justy <- gtable_filter(ggplotGrob(justy), 'axis-l|ylab', trim=F)
justy <- ggdraw(justy)

## permutation plot
perm <- ggplot(ld.sig, aes(x=iter, y=coef)) +
            geom_hline(yintercept=0, lty=2, lwd=1, colour="grey60") +
            geom_errorbar(aes(ymin=lwr, ymax=upr), lwd=0.5, colour="#8da0cb", width=0) +
            geom_errorbar(aes(ymin=coef - cse, ymax=coef + cse), lwd=0.5, colour="#fc8d62", width=0) +
            geom_point(size=0.5, pch=21, fill="red4", color="red4") +
            ylab("+1\u0394\u2103 on \u0394 Mental Health Issues\n(Percentage Points)") +
            xlab("Long-Diff. Model Permutation") +
            coord_cartesian(ylim=c(-2.5,5)) +
            scale_x_continuous(breaks = c(0,20,40,60)) +
            theme(plot.title = element_text(face="bold", size=40),
                  axis.title.x = element_text(size=30, angle=0, face="bold"),
                  axis.title.y = element_blank(),
                  axis.text.y = element_blank() ,
                  axis.text.x = element_text(size=25, face="bold"),
                  panel.grid.major = element_blank(),
                  panel.grid.minor = element_blank(),
                  panel.background = element_blank(),
                  axis.line.x = element_line(),
                  axis.line.y = element_blank(),
                  axis.ticks.y = element_blank(),
                  legend.title = element_blank()
            )
hist <- axis_canvas(perm, axis = "y", coord_flip=TRUE) + 
        geom_histogram(data = ldcomb, aes(x = hist.coef), bins=20, fill="#8da0cb", color='gray60', alpha=0.5, size=0.1) +
        coord_flip(xlim=c(-2.5,5)*-1) + # need to multiply by negative one to keep it consistent
        scale_x_reverse()
perm <- insert_yaxis_grob(perm, hist, position = "right", width = grid::unit(0.1, "null"))
perm <- ggdraw(perm)

# plot----
cairo_pdf(filename = './figure3.pdf', width=16, height=8, fallback_resolution = 1200)
grid.arrange(arrangeGrob(map, p, ncol=2), 
             arrangeGrob(rectGrob(gp=gpar(col="white")), legend1, rectGrob(gp=gpar(col="white")), legend2, ncol=4, widths=c(0.05, 0.45, 0.05, 0.45)),
             arrangeGrob(rectGrob(gp=gpar(col="white")), 
                         arrangeGrob(justy, coef, perm, ncol=3, widths=c(0.08, 0.44, 0.44)), 
                         rectGrob(gp=gpar(col="white")), ncol=3, widths=c(0.02, 0.96, 0.02)),
             nrow=3, 
             heights=c(0.52, 0.08, 0.4)) # plot all together
grid.text("a", x=unit(0.02, "npc"), y=unit(0.98, "npc"), gp=gpar(fontsize=40,font=2))
grid.text("b", x=unit(0.53, "npc"), y=unit(0.98, "npc"), gp=gpar(fontsize=40,font=2))
grid.text("c", x=unit(0.13, "npc"), y=unit(0.15, "npc"), gp=gpar(fontsize=40,font=2))
grid.text("d", x=unit(0.56, "npc"), y=unit(0.15, "npc"), gp=gpar(fontsize=40,font=2))
dev.off()
```


```{r fig.longdiff, out.width="100%", fig.cap="Longer-term warming and rates of mental health issues. Panel (a) presents the change in prevalence of mental health issues between 2002-2006 and 2007-2011. Points are scaled by the number of respondents from each city. Panel (b) plots the change in grid cell average maximum temperatures over this same period. Panel (c) plots the coefficient estimates from the long-differences model including all observations, the model excluding cities missing any years, and the model including only the cities missing one or more years in our sample, respectively. Panel (d) plots the results of a permutation of all possible long-differences combinations in our data where the difference between periods is at least five years. The coefficient estimates from this process whose 95\\% confidence intervals do not contain zero are points in panel (d). We display all coefficient estimates from this permutation process in the histogram on the right of panel (d). Orange error lines represent one standard error; blue lines represent 95\\% confidence intervals."}
knitr::include_graphics("figure3.pdf")
```

```{r}
longdiffs <- fig_num(name = "longdiffs", caption = "")
```

# Figure 4

```{r fig4, eval=FALSE}
# spring----
ld <- data.table(coef = c(ld.spr$coef[2], ld.full.spr$coef[2], ld.partial.spr$coef[2]), se = c(ld.spr$cse[2], ld.full.spr$cse[2], ld.partial.spr$cse[2]), model = c("All", "Full", "Missing"))
ld[, upr := coef + 1.96*se]
ld[, lwr := coef - 1.96*se]
coef.spr <- ggplot(ld, aes(model, coef)) +
        geom_hline(yintercept=0, lty=2, lwd=1, colour="grey60") +
        geom_errorbar(aes(ymin=lwr, ymax=upr), lwd=2, colour="#8da0cb", width=0) +
        geom_errorbar(aes(ymin=coef - se, ymax=coef + se), lwd=5, colour="#fc8d62", width=0) +
        geom_point(size=8, pch=21, fill="red4") +
        ylab("+1\u0394\u2103 on \u0394 Mental Health Issues\n(Percentage Points)") +
        xlab("Spring") +
        coord_cartesian(ylim=c(-4,4)) +
        theme(plot.title = element_text(face="bold", size=40),
              axis.title.x = element_text(size=30, angle=0, face="bold"),
              axis.title.y = element_blank(),
              axis.text.y = element_blank(),
              axis.text.x = element_text(size=25, face="bold"),
              panel.grid.major = element_blank(),
              panel.grid.minor = element_blank(),
              panel.background = element_blank(),
              axis.line.x = element_line(),
              axis.line.y = element_blank(),
              axis.ticks.y = element_blank(),
              legend.title = element_blank()
        )

# summer----
ld <- data.table(coef = c(ld.sum$coef[2], ld.full.sum$coef[2], ld.partial.sum$coef[2]), se = c(ld.sum$cse[2], ld.full.sum$cse[2], ld.partial.sum$cse[2]), model = c("All", "Full", "Missing"))
ld[, upr := coef + 1.96*se]
ld[, lwr := coef - 1.96*se]
coef.sum <- ggplot(ld, aes(model, coef)) +
        geom_hline(yintercept=0, lty=2, lwd=1, colour="grey60") +
        geom_errorbar(aes(ymin=lwr, ymax=upr), lwd=2, colour="#8da0cb", width=0) +
        geom_errorbar(aes(ymin=coef - se, ymax=coef + se), lwd=5, colour="#fc8d62", width=0) +
        geom_point(size=8, pch=21, fill="red4") +
        ylab("+1\u0394\u2103 on \u0394 Mental Health Issues\n(Percentage Points)") +
        xlab("Summer") +
        coord_cartesian(ylim=c(-4,4)) +
        theme(plot.title = element_text(face="bold", size=40),
              axis.title.x = element_text(size=30, angle=0, face="bold"),
              axis.title.y = element_blank(),
              axis.text.y = element_blank(),
              axis.text.x = element_text(size=25, face="bold"),
              panel.grid.major = element_blank(),
              panel.grid.minor = element_blank(),
              panel.background = element_blank(),
              axis.line.x = element_line(),
              axis.line.y = element_blank(),
              axis.ticks.y = element_blank(),
              legend.title = element_blank()
        )


# fall----
ld <- data.table(coef = c(ld.fall$coef[2], ld.full.fall$coef[2], ld.partial.fall$coef[2]), se = c(ld.fall$cse[2], ld.full.fall$cse[2], ld.partial.fall$cse[2]), model = c("All", "Full", "Missing"))
ld[, upr := coef + 1.96*se]
ld[, lwr := coef - 1.96*se]
coef.fall <- ggplot(ld, aes(model, coef)) +
        geom_hline(yintercept=0, lty=2, lwd=1, colour="grey60") +
        geom_errorbar(aes(ymin=lwr, ymax=upr), lwd=2, colour="#8da0cb", width=0) +
        geom_errorbar(aes(ymin=coef - se, ymax=coef + se), lwd=5, colour="#fc8d62", width=0) +
        geom_point(size=8, pch=21, fill="red4") +
        ylab("+1\u0394\u2103 on \u0394 Mental Health Issues\n(Percentage Points)") +
        xlab("Fall") +
        coord_cartesian(ylim=c(-4,4)) +
        theme(plot.title = element_text(face="bold", size=40),
              axis.title.x = element_text(size=30, angle=0, face="bold"),
              axis.title.y = element_blank(),
              axis.text.y = element_blank(),
              axis.text.x = element_text(size=25, face="bold"),
              panel.grid.major = element_blank(),
              panel.grid.minor = element_blank(),
              panel.background = element_blank(),
              axis.line.x = element_line(),
              axis.line.y = element_blank(),
              axis.ticks.y = element_blank(),
              legend.title = element_blank()
        )

# winter----
ld <- data.table(coef = c(ld.win$coef[2], ld.full.win$coef[2], ld.partial.win$coef[2]), se = c(ld.win$cse[2], ld.full.win$cse[2], ld.partial.win$cse[2]), model = c("All", "Full", "Missing"))
ld[, upr := coef + 1.96*se]
ld[, lwr := coef - 1.96*se]
coef.win <- ggplot(ld, aes(model, coef)) +
        geom_hline(yintercept=0, lty=2, lwd=1, colour="grey60") +
        geom_errorbar(aes(ymin=lwr, ymax=upr), lwd=2, colour="#8da0cb", width=0) +
        geom_errorbar(aes(ymin=coef - se, ymax=coef + se), lwd=5, colour="#fc8d62", width=0) +
        geom_point(size=8, pch=21, fill="red4") +
        ylab("+1\u0394\u2103 on \u0394 Mental Health Issues\n(Percentage Points)") +
        xlab("Winter") +
        coord_cartesian(ylim=c(-4,4)) +
        theme(plot.title = element_text(face="bold", size=40),
              axis.title.x = element_text(size=30, angle=0, face="bold"),
              axis.title.y = element_blank(),
              axis.text.y = element_blank(),
              axis.text.x = element_text(size=25, face="bold"),
              panel.grid.major = element_blank(),
              panel.grid.minor = element_blank(),
              panel.background = element_blank(),
              axis.line.x = element_line(),
              axis.line.y = element_blank(),
              axis.ticks.y = element_blank(),
              legend.title = element_blank()
        )

# justy----
justy <- ggplot(ld, aes(model, coef)) +
        geom_hline(yintercept=0, lty=2, lwd=1, colour="grey60") +
        geom_errorbar(aes(ymin=lwr, ymax=upr), lwd=2, colour="#8da0cb", width=0) +
        geom_errorbar(aes(ymin=coef - se, ymax=coef + se), lwd=5, colour="#fc8d62", width=0) +
        geom_point(size=8, pch=21, fill="red4") +
        ylab("Effect of +1\u2103 on \n\u0394 Mental Health Issues\n(Percentage Points)") +
        xlab("Long-Diff. Model") +
        coord_cartesian(ylim=c(-4,4)) +
        theme(plot.title = element_text(face="bold", size=40),
              axis.title.x = element_text(size=30, angle=0, face="bold"),
              axis.title.y = element_text(size=35, angle=90, vjust=-0.25, face="bold"),
              axis.text.y = element_text(size=25, angle=0, face="bold") ,
              axis.text.x = element_text(size=25, face="bold"),
              panel.grid.major = element_blank(),
              panel.grid.minor = element_blank(),
              panel.background = element_blank(),
              axis.line.x = element_line(),
              axis.line.y = element_line(),
              legend.title = element_blank()
        )
justy <- gtable_filter(ggplotGrob(justy), 'axis-l|ylab', trim=F)
justy <- ggdraw(justy)

# plot----
cairo_pdf(filename = './figure4.pdf', width=16, height=6.5, fallback_resolution = 1200)
grid.arrange(justy, arrangeGrob(coef.spr, coef.sum, coef.fall, coef.win, ncol=4), ncol=2, widths=c(0.12, 0.88))
grid.text("a", x=unit(0.14, "npc"), y=unit(0.96, "npc"), gp=gpar(fontsize=45,font=2))
grid.text("b", x=unit(0.36, "npc"), y=unit(0.96, "npc"), gp=gpar(fontsize=45,font=2))
grid.text("c", x=unit(0.58, "npc"), y=unit(0.96, "npc"), gp=gpar(fontsize=45,font=2))
grid.text("d", x=unit(0.8, "npc"), y=unit(0.96, "npc"), gp=gpar(fontsize=45,font=2))
dev.off()
```


```{r fig.ldseasons, out.width="100%", fig.cap="Multiyear warming and prevalence of mental health problems across the seasons of the year. Panel (a) plots the results of estimating the long-differences model for spring, panel (b) for summer, panel (c) for fall, and panel (d) for winter. Orange error lines represent one standard error; blue lines represent 95\\% confidence intervals."}
knitr::include_graphics("figure4.pdf")
```

```{r}
seasons <- fig_num(name = "seasons", caption = "")
```

# Figure 5

```{r fig5, eval=FALSE}
# set up diff-n-diff variables----
# all cities that fell in FEMA disaster zones for Katrina
brfss[, treat := (state %like% "LA|MS") | (name %like% "Birmingham AL|Mobile AL|Tuscaloosa AL|Miami FL|Key West FL|Naples FL|Panama City FL|Pensacola FL")] 

# set pre-Katrina as pre and post-Katrina as post
brfss[, post := time > "2005-08-24"]

# run randomization inference over 100K samples----
# set.seed(1987) # set seed for replicable plot.
# bootstrap <- foreach (i = 1:100000) %dopar% {
#                 brfss[, treat.fake := sample(treat, size=length(treat), replace=TRUE)]
#                 fake <- felm(any_ment ~  post*treat.fake
#                              | 0
#                              | 0
#                              | state, data=brfss, exactDOF=FALSE)
#                 coefs <- as.data.table(cbind(i,fake$coef[4],fake$cse[4]))
#                 setnames(coefs, c('iteration','coef','cse'))
#                 print(paste(i))
#                 return(coefs)
#             }
# placebo.coefs <- rbindlist(bootstrap)
# saveRDS(placebo.coefs, "./randomization_inference_coefficients.rds", compress=TRUE)
placebo.coefs <- readRDS("./randomization_inference_coefficients.rds")

# calculate estimate of true effect----
true <- felm(any_ment ~ post*treat | 0 | 0 | state, data=brfss, exactDOF=TRUE)
true <- as.data.table(cbind(true$coef[4], true$cse[4]))
setnames(true, c('coef','cse'))

# plots----
## permutation plot
p1 <- ggplot() +
        geom_histogram(data=placebo.coefs, aes(x=coef), fill="gray70", size=0.25, alpha=0.75, bins=100) +
        geom_vline(xintercept = true$coef, color="#fc8d62", linetype='dashed', size=1.5) +
        geom_text(aes(x = true$coef-1.3, y=700, label="Katrina\nCoefficient\u2192"), size=8.5) +
        geom_text(aes(x = median(placebo.coefs$coef), y=700, label="Permutation\nCoefficients"), size=8.5) +
        xlab('\u0394 P(Mental Health Issue), Diff. in Diff.') +
        ylab('Frequency') +
        coord_cartesian(xlim=c(-2.5,5.5), expand=FALSE) +
        theme(title = element_text(size=20, angle=0, vjust=.5, face="bold"),
              plot.title = element_text(hjust=0.025),
              axis.title.x = element_text(size=22, angle=0, vjust=-0, face="bold"),
              axis.title.y = element_text(size=25,angle=90, vjust=1, face="bold"),
              axis.text.y = element_text(size=20,angle=0, vjust=.5, face="bold") ,
              axis.text.x = element_text(size=20, face="bold"),
              panel.grid.major = element_blank(),
              panel.grid.minor = element_blank(),
              panel.background = element_blank(),
              axis.line.x = element_line(),
              axis.line.y = element_line(),
              legend.title = element_blank(),
              legend.position = c(.5, .15),
              legend.text = element_text(size=25),
              legend.key = element_rect(colour = "black", fill="transparent"),
              legend.box = "vertical"
        )

## diff-n-diff mean plot
sem <- function(x) {sqrt(var(x, na.rm = TRUE) / length(na.omit(x)))} # sem function
diff <- brfss[!is.na(post), .(any_ment = mean(any_ment, na.rm=T), any_ment_se = sem(any_ment)), by=.(post, treat)]
diff[, upr := any_ment + any_ment_se]
diff[, lwr := any_ment - any_ment_se]
diff[, treat := as.character(treat)]
diff[, post := as.character(post)]
diff[treat=="TRUE" , treat := "Katrina-Exposed"]
diff[treat=="FALSE", treat := "Not-Katrina-Exposed"]
diff[post=="TRUE", post := "Post-Period"]
diff[post=="FALSE", post := "Pre-Period"]
diff[, post := as.factor(post)]
diff[, post := relevel(post, ref="Pre-Period")]
diff[, treat := as.factor(treat)]

dodge <- position_dodge(width=0.90)
p2 <- ggplot(data = diff, aes(y = any_ment, x = factor(treat), fill=factor(post))) +
          geom_bar(stat = "identity", position=dodge, color='black') +
          geom_linerange(aes(ymax = upr, ymin = lwr), position=dodge) +
          scale_fill_manual(values=c("#8da0cb","#fc8d62"), breaks=c("Pre-Period","Post-Period"), labels=c("Pre-Period","Post-Period")) +
          coord_cartesian(xlim = c(0.5, 2.5), ylim = c(0,42), expand=FALSE) +
          ylab("% Mental Health Issues") +
          xlab("Treatment Group") +
          geom_hline(yintercept = 0) +
          ggplot2::annotate("text",x=1,y=35.75, size=9, label="***") +
          ggplot2::annotate("text",x=2,y=35.75, size=9, label="*") +
          ggplot2::annotate("text",x=1.5,y=38.75, size=9, label="***") +
          geom_segment(aes(x = 0.75, y = 35.5, xend = 1.25, yend = 35.5), size=1) +
          geom_segment(aes(x = 0.75, y = 35.25, xend = 0.75, yend = 35.5), size=1) +
          geom_segment(aes(x = 1.25, y = 35.25, xend = 1.25, yend = 35.5), size=1) +
          geom_segment(aes(x = 1.75, y = 35.5, xend = 2.25, yend = 35.5), size=1) +
          geom_segment(aes(x = 1.75, y = 35.25, xend = 1.75, yend = 35.5), size=1) +
          geom_segment(aes(x = 2.25, y = 35.25, xend = 2.25, yend = 35.5), size=1) +
          geom_segment(aes(x = 1, y = 38.5, xend = 2, yend = 38.5), size=1) +
          geom_segment(aes(x = 1, y = 38.25, xend = 1, yend = 38.5), size=1) +
          geom_segment(aes(x = 2, y = 38.25, xend = 2, yend = 38.5), size=1) +
          theme(axis.title.x = element_text(size=22, angle=0, vjust=-0, face="bold"),
                axis.title.y = element_text(size=25,angle=90, vjust=1, face="bold"),
                axis.text.y = element_text(size=20,angle=0, vjust=.5, face="bold") ,
                axis.text.x = element_text(size=20, face="bold"),
                panel.grid.major = element_blank(),
                panel.grid.minor = element_blank(),
                panel.background = element_blank(),
                axis.line = element_line(),
                legend.position=c(0.06,0.15),
                legend.background = element_rect("white"),
                legend.title = element_blank(),
                legend.text=element_text(size=25, face="bold"),
                legend.key = element_rect(colour = "black", fill="white")
                )

# plot----
cairo_pdf(filename = "./figure5.pdf", width=14, height=4.6, fallback_resolution = 1200)
grid.arrange(p2, p1, ncol=2)
grid.text("a", x=unit(0.08, "npc"), y=unit(0.96, "npc"), gp=gpar(fontsize=40,font=2))
grid.text("b", x=unit(0.605, "npc"), y=unit(0.96, "npc"), gp=gpar(fontsize=40,font=2))
dev.off()
```


```{r fig.katrina, out.width="100%", fig.cap="Worsened mental health in areas affected by Hurricane Katrina. Panel (a) depicts the mean prevalence of mental health difficulties for the Katrina-exposed and non-Katrina-exposed groups before and after Hurricane Katrina landfall. Locales not exposed to Katrina observed a drop in rates of mental health issues, while locales exposed to Katrina experienced an increase. The difference-in-differences estimate indicates Hurricane Katrina exposure associates with an increase in prevalence of mental health issues. All differences are significant at $\\alpha = 0.05$. Error bars are SEM. Panel (b) plots the results of randomly permuting 100,000 treatment status assignments and conducting our difference-in-differences regression for each permutation. The true difference-in-differences estimate falls substantially outside the distribution of permutation coefficients."}
knitr::include_graphics("figure5.pdf")
```

```{r}
katrina <- fig_num(name = "katrina", caption = "")
```
